// Mood Analysis
// BY Jared Bauer, Alex Jansen and Jesse Cirimele
// Based on example code by Krister Olsson <http://www.tree-axis.com> 

// Displays realtime mood  
// and performing real-time spectrum analysis. 
// Click to start/stop

// Created 1 May 2005
// Updated 20 August 2011

import krister.Ess.*;
import java.util.Arrays;
import java.io.*;

import org.apache.commons.math.transform.FastFourierTransformer;
import org.apache.commons.math.analysis.*;
import org.apache.commons.math.complex.*;


AudioInput myChannel;
//old fft library
krister.Ess.FFT myFFT;
//new fft library
FastFourierTransformer transformer;

PFont font;

int bufferSize;
int bufferDuration;
boolean playing = false;
boolean meanFlag = false;
float mean = 0;
PrintWriter output;
int startTime;

int timeCounter = 1;

//this sets the looping framerate of the drawing portion of the application
int OUR_FRAMERATE = 32;
//this sets how often we log the output of the program (in seconds)
int OUR_LOG_RATE = 30;
//this is where our logged data and program output is saved
String DATA_LOG_FILE = "data/tension-energy-"+year()+month()+day()+"_"+hour()+minute()+second()+".csv";

//this variable sets the volume threshold to detect speech.
float volume_threshold = .015;
//these variables are globals used to track history for averaging (not actually 5 seconds)
float[] tension_history;
float[] energy_history;
//this variable sets how long our window to average over is
int HISTORY_SIZE = OUR_FRAMERATE*OUR_LOG_RATE;
//this variable tracks where we are in the array
int current_history_location;

int SAMPLE_RATE = 44100;

void setup() {
  size(540,220);
  startTime = millis();
  
  font = loadFont("CourierNew36.vlw"); 
  
  //create a logger file for the output data
  output = createWriter(DATA_LOG_FILE);
  
  // start up Ess
  Ess.start(this);
  int sampleRate = SAMPLE_RATE;
  bufferSize=2048; 
 
  myChannel=new AudioInput(bufferSize);
  bufferDuration=myChannel.ms(sampleRate*2);
  
   myChannel.sampleRate(sampleRate);

  //old FFT library
  myFFT=new krister.Ess.FFT(1024);
  //new FFT library
  transformer = new FastFourierTransformer();

  frameRate(OUR_FRAMERATE);
  noSmooth();

  //setup history for rms (volume) and pitch
  tension_history = new float[HISTORY_SIZE];
  energy_history = new float[HISTORY_SIZE];
  for(int i = 0; i < HISTORY_SIZE; i++) {
    tension_history[i] = 0;
    energy_history[i] = 0;
  }
  //reset the current place in array to 0
  current_history_location = 0;
}

void draw() {
  background(25);
  drawSpectrum();
  drawSamples();
  drawLocation();
  
}

void drawSpectrum() {
  noStroke();

  //trying to test new fft code
  double theBuffer[];
  Complex theBuffer2[];
  theBuffer = convertFloatsToDoubles( myChannel.buffer2 );
  theBuffer2 = convertFloatsToComplexs( myChannel.buffer2 );
  
  //try low-pass filter
  float[] signal_lp = lowPassFilter( myChannel.buffer2, 500 );
  
  //Complex result[] = transformer.transform(theBuffer);
  Complex[] theBuffer3 = convertFloatsToComplexs( signal_lp );
  Complex result[] = FFT.fft( theBuffer3 );

  for (int i=0; i<256; i++) {
    double magnitudeFFT = sqrt((float)(result[i].re()*result[i].re() + result[i].im()*result[i].im()));
    float temp=max(0,185-((int)magnitudeFFT));

    rect(i,temp+.5,1,height-temp+.5);
    stroke(255,153,0);
  }
}

void drawSamples() {
  stroke(255);
  
  //try low-pass filter
  float[] signal_lp = lowPassFilter( myChannel.buffer2, 500 );

  // interpolate between 0 and writeSamplesSize over writeUpdateTime
  int interp=(int)max(0,(((millis()-myChannel.bufferStartTime)/(float)bufferDuration)*bufferSize));

  for (int i=0;i<256;i++) {
    float left=100;
    float right=100;

    //plot full signal
    /*if (i+interp+1<myChannel.buffer2.length) {
      left-=myChannel.buffer2[i+interp]*75.0;
      right-=myChannel.buffer2[i+1+interp]*75.0;
    }*/
    
    //plot low pass signal
    if (i+interp+1<signal_lp.length) {
      left-=signal_lp[i+interp]*75.0;
      right-=signal_lp[i+1+interp]*75.0;
    }

    line(i,left,i+1,right);
    stroke(255,153,0);  
  }
}

//Draw a grid that a dot will appear on to indicate where the audio is in our acquistic space
void drawLocation() {
    int xOrigin = 425;
    int yOrigin = 110;
    float sumOfSpectrums = 0;
    float sumOfSpectrumsSq = 0;
    float energy = 0;
    float ampVariance = 0;
    float maxFreq = 0;
    float maxBand = 0;
    int maxBandLoc =  0;
    int m = millis();
    int[] maxBandArray = new int[100];
    int sum = 0;    
    
    fill(45);
    rect(330, 5, 190, 190);
    
    stroke(255,153,0);
    noStroke();
    
    textFont(font, 14); 
    fill(255,153,0);
    textAlign(RIGHT);
    text("Energetic", 328, 30); 
    text("Tired", 328, 192);
    textAlign(CENTER);
    text("Calm", 351, 210);
    text("Tense", 490, 210);
  
  int spectrumLength = myFFT.spectrum.length;

  //PITCH DETECTION
  float thePitch = detectedPitchFrequency( myChannel.buffer2 );
 
 //CALCULATE RMS for intensity
  float rms_sum = 0;
  for(int i = 0; i < myChannel.buffer2.length; i++) {
     rms_sum += myChannel.buffer2[i]*myChannel.buffer2[i];
  }
  float rms = sqrt(rms_sum / myChannel.buffer2.length);
  
  //SETUP: NEW GRAPHING
  //these weighting factors are used to scale the effects of volume and pitch
  float weighting_rms_energy = .7275;
  float weighting_rms_tension = .605;
  float weighting_pitch_energy = .3075;
  
  //these averages are used to determine loud vs soft and high vs low for participant speech
  float rms_baseline = .06; //this is roughly the average volume we expect
  float pitch_baseline = 168; //this is roughly the average pitch we expect
  float min_pitch = 72; //this is roughly the minimum pitch we expect
  
  //get the current energy and tension to put into the history
  float energy_current = 0;
  float tension_current = 0;
  if (rms > volume_threshold) {
    //volume: rms_energy + pitch_energy
    float pitch_energy = ((thePitch - pitch_baseline)/(pitch_baseline-min_pitch)) * weighting_pitch_energy;
    float rms_energy = ((rms - rms_baseline)/(rms_baseline-volume_threshold)) * weighting_rms_energy;
    energy_current = rms_energy + pitch_energy;
    if(energy_current > 1) energy_current = 1;
    if(energy_current < -1) energy_current = -1;

    //tension: rms intensity
    tension_current = rms_energy * weighting_rms_tension;
    if(tension_current > 1) tension_current = 1;
    if(tension_current < -1) tension_current = -1;
  }
  
  //set our current values so that we can average
  tension_history[current_history_location] = tension_current;
  energy_history[current_history_location] = energy_current;
  
  float tension_average_x = get_avg_history(tension_history);
  float energy_average_y = get_avg_history(energy_history);
  
  current_history_location++;
  if (current_history_location >= HISTORY_SIZE)
    current_history_location = 0;
  
  float tension_plot_x = map ( tension_average_x, -1.1, 1.1, -80, 80 );
  float energy_plot_y = map ( energy_average_y, -1.1, 1.1, -80, 80 );
  
   //LOGGING
   //write to the file every 3 seconds
  if ((m - startTime) / 3000 > timeCounter ){
    println(DATA_LOG_FILE);
//    output = createWriter(DATA_LOG_FILE);
    output.println("timestamp, " + hour()+":"+minute()+":"+second()+ ", rms , " + rms + ", pitch, " + thePitch + ", Tension_AVG, " + tension_average_x + ", Energy_AVG, " + energy_average_y);
    println( " Tension , " + tension_average_x + ", Energy, " + energy_average_y);
    println( " rms " + rms + " freq " + thePitch + " -- tension " + tension_current + " energy " + energy_current + " -- time location " +   current_history_location);
    
    output.flush();
    
    timeCounter++;
  }
  ellipse(xOrigin + tension_plot_x*1.8, yOrigin - energy_plot_y*1.8, 10,10);  
}

float get_avg_history(float[] history_array) {
  float temp = 0;
  int the_length = history_array.length;
  for (int i = 0; i < the_length; i++) {
    temp += history_array[i];
  }
  temp = temp / (float)the_length;
  return temp;
}

void mousePressed(){
  
 if (playing == false) {
    myChannel.start();
    playing = true;
    //exit;
 }
 else if (playing == true) {
    myChannel.stop();
    playing = false;
    output.close();
    exit();
 }
}
// clean up Ess before exiting

public void stop() {
  Ess.stop();
  super.stop();
}

//helper code to convert float[] to double[]
public static double[] convertFloatsToDoubles(float[] input)
{
    if (input == null)
    {
        return null; // Or throw an exception - your choice
    }
    double[] output = new double[input.length];
    for (int i = 0; i < input.length; i++)
    {
        output[i] = input[i];
    }
    return output;
}

//helper code to convert float[] to Complex[]
public static Complex[] convertFloatsToComplexs(float[] input)
{
    if (input == null)
    {
        return null; // Or throw an exception - your choice
    }
    Complex[] output = new Complex[input.length];
    for (int i = 0; i < input.length; i++)
    {
        output[i] = new Complex(input[i],0);
    }
    return output;
}

//helper code to convert Complex[] to float[]
public static float[] convertComplexsToFloats(Complex[] input)
{
    if (input == null)
    {
        return null; // Or throw an exception - your choice
    }
    float[] output = new float[input.length];
    for (int i = 0; i < input.length; i++)
    {
        output[i] = (float) input[i].re();
    }
    return output;
}

// this is the autocorrelation-like function. it compares a signal to itself shifted by
// amounts from 0 to M_max. N (the signal length) is the length of sounds in question.
// 
public float[] AMDFunction ( float[] signal) {
  //this sets the maximum shift value. at 44100 sample rate for 100Hz it should be around 450;
  float lowest_frequency = 80; //lowest frequency in Hz
  int M_max = (int)((1/lowest_frequency)*SAMPLE_RATE); //number of samples
  
  //set first 20 samples
  for( int i=0; i<20; i++ ) {
    signal[i] = 0;
  }
  
  int N = signal.length;
  if(M_max > N)
    M_max = N;
  float[] AMDF_arr = new float[M_max];
  
  for (int m = 20; m < M_max; m++) {
    AMDF_arr[m] = 0.0;
    for (int i=0; i<N; i++ ) {
      int i_m = i - m;
      //we loop the array for values outside the array
      if( (i_m) < 0 ) {
        i_m += N;
      }
      AMDF_arr[m] += abs( signal[i] - signal[i_m] );
    }
    AMDF_arr[m] = AMDF_arr[m] / N;
  }
  
  return AMDF_arr;
}

public void AMDF_print (double[] AMDF_arr) {
  int N = AMDF_arr.length;
  for (int i=0; i < N; i++) {
    print(" "+AMDF_arr[i]);
  }
  println(" ");
}

// this function lets all frequencies lower than frequency_lp through
// but frequencies higher are removed. cut_frequency should be in Hz.
public float[] lowPassFilter (float[] signal, float frequency_lp) {
  Complex[] signal_complex = convertFloatsToComplexs( signal );
  //convert to the frequency domain
  Complex[] freq_domain = FFT.fft( signal_complex );
  
  // generate a low-pass filter of the spectrum
  int SIZE = freq_domain.length;
  Complex[] filterSpectrum = new Complex[SIZE];
  for (int i = 0; i < SIZE; i++)
  {
      if (i < (SIZE * frequency_lp / SAMPLE_RATE))
      {
          // passband
          filterSpectrum[i] = new Complex(freq_domain[i].re()*1.0, freq_domain[i].im()*0.0);
      }
      else
      {
          // stopband
          filterSpectrum[i] = new Complex(0.0,0.0);
      }
  }
  
  //convert back time-domain
  Complex[] filteredSignal = FFT.ifft( filterSpectrum );
  
  //convert back to floats from Complexs
  float[] filteredSignal_float = convertComplexsToFloats( filteredSignal );
  
  //return final signal array
  return filteredSignal_float;
}  

public float detectedPitchFrequency ( float[] signal ) {
  //apply a low pass filter (human voice is between 100 and 900 Hz)
  float [] signal_lp = lowPassFilter ( signal, 900 );
  
  // Center clipping helps make the pitch detection more accurate. this should
  // retuce octive errors
  //not doing this yet
  
  //AMDF (Average Magnitude Difference Function)
  float[] AMDF_arr = AMDFunction( signal_lp );
    
  //Pitch = Index of min (AMDF)
  int index = getMin( AMDF_arr, 20 );
  
  //convert index to frequency. index is in samples.
  float cycleTime = (float)index/(float)SAMPLE_RATE;
  float finalFrequency = 1.0/cycleTime;
  
  return finalFrequency;
}

public int getMin( float[] theArray, int startIndex ) {
  int N = theArray.length;
  int minIndex = 0;
  float theMin = 50000;
  for (int i= startIndex; i < N; i++) {
    if( theArray[i] < theMin ) {
      theMin = theArray[i];
      minIndex = i;
    }
  }
  return minIndex;
}
